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We develop and analyze a two-mode phase-field-crystal model to describe fee ordering. The 
model is formulated by coupling two different sets of crystal density waves corresponding to (111) 
and (200) reciprocal lattice vectors, which are chosen to form triads so as to produce a simple free- 
energy landscape with coexistence of crystal and liquid phases. The feasibility of the approach is 
demonstrated with numerical examples of polycrystalline and (111) twin growth. We use a two-mode 
amplitude expansion to characterize analytically the free-energy landscape of the model, identifying 
parameter ranges where fee is stable or metastable with respect to bee. In addition, we derive 
analytical expressions for the elastic constants for both fee and bee. Those expressions show that a 
non-vanishing amplitude of [200] density waves is essential to obtain mechanically stable fee crystals 
with a non-vanishing tetragonal shear modulus (Cn — Ci2)/2. We determine the model parameters 
for specific materials by fitting the peak liquid structure factor properties and solid density wave 
amplitudes following the approach developed for bcc [K.-A. Wu and A. Karma, Phys. Rev. B 76, 
184107 (2007)]. This procedure yields reasonable predictions of elastic constants for both bcc Fe 
and fee Ni using input parameters from molecular dynamics simulations. The application of the 
model to two-dimensional square lattices is also briefiy examined. 



PACS numbers: 61.72.Mm,68.08.De,68.08.-p,81.16.Rf 

I. INTRODUCTION AND SUMMARY 

The phase-field-crystal (PFC) method has emerged as 
an attractive computational approach to simulate the 
evolution of crystalline patterns [TH^. By resolving the 
crystal density field, it naturally incorporates defects 
and elastic interactions arising from localized and large 
scale distortions of this field, respectively. Moreover, this 
method can in principle be used to simulate microstruc- 
tural evolution on diffusive time scales that are much 
longer than typical time scales accessible by molecular 
dynamics (MD) simulations. 

Like classical density function theory (DFT) , the PFC 
method is based on representing the free-energy of a ma- 
terial by a functional of its density [7l-fT3]. However, clas- 
sical DFT and the PFC method use different functionals 
to achieve different goals. Classical DFT seeks a phys- 
ically realistic mean-field description of the crystal den- 
sity field n{r) to reproduce quantitatively as accurately 
as possible the properties of a material. Since n{r) is 
sharply peaked around mean atomic positions, this gen- 
erally requires a very large number of terms in the tra- 
ditional expansion of the number density as a sum of 
density waves 




where each Ki represents a different reciprocal lattice 
vector (RLV) in this unrestricted sum. In contrast, by us- 
ing a considerably simplified density functional, the PFC 
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method essentially restricts this sum to a much smaller 
set of reciprocal lattice vectors in order to simulate effi- 
ciently the evolution of the crystal field on length and 
time scales as large as possible. Recent studies have 
shown that, despite this loss of realism, PFC models 
are able to reproduce quantitatively certain key proper- 
ties that influence microstructural evolution such as the 
crystal-melt interfacial free-energy [HI US], the bulk mod- 
ulus [15], and grain-boundary energies [TS], which have 
been computed for the test case of pure Fe. 

Despite this progress, the PFC method has only been 
developed for a small set of crystal structures. The orig- 
inal formulation of Elder et al. pQ |2] uses the same free- 
energy functional as the Swift-Hohenberg model of pat- 
tern formation [161 117) of the form 




with the free-energy density 

f=t[a + X{q^+^^f]^ + g^, (3) 

where <j) represents the crystal density field. This one- 
mode model essentially truncates the sum ([!]) to one set 
of RLVs with equal magnitude \Ki\ = go since higher 
K modes have much smaller amplitude. As a result, it 
favors crystal structures for which the principal RLVs 
can form "triads" (i.e. closed triangles), which include 
hexagonal and body-centered-cubic (bcc) ordering in two 
and three dimensions, respectively. Aside from favoring 
those structures, triad interactions are essential for solid- 
liquid coexistence. This is because in a weakly nonlinear 
expansion of the bulk free-energy density of the form, 
/ = C2U^ + cau'^ -|- -f . . . (with Ui — u for all principal 
RLVs), triads contribute a cubic term with a negative 
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coefficient C3 < 0. Since C2 and C4 are both positive, this 
cubic term is responsible for the existence of a free-energy 
barrier between the two minima of / corresponding to 
hquid {u = 0) and sohd {us > 0). 

In this paper, we use a "two-mode" phase-field-crystal 
model to model face-centered-cubic (fee) structures, 
which has the free-energy density 

/ = ^ [a + A(V^ + qiraw' + qlf + n)] + g^. (4) 

This model truncates the sum ([ij to two sets of RLVs 
with magnitude \Ki\ — go and \K'^\ = qi, respectively, 
where the first set corresponds in general to the princi- 
pal RLVs and the second to some other set with larger 
wavevector magnitude; all other RLVs have much smaller 
amplitude. This construct provides more flexibility to 
form triad interactions by combining RLVs from those 
two sets, and hence to describe other crystal structures. 
We demonstrate this here for fee ordering, which is ob- 
tained by choosing the sets {Ki} and {K'^} to correspond 
to (111) (principal set) and (200) RLVs, respectively, 
with qi/qo = \/4/3. While there is in principle freedom 
in the choice of the second set {i?^} for a given structure, 
we have chosen this set such that qi > qo is as small as 
possible, as desired for computational efficiency. 

The form Q reduces in the limit ri = to the 
free-energy density introduced by Lifshitz and Petrich 
[H] as a generalization of the Swift-Hohenberg model 
to describe two-dimensional quasiperiodic patterns ob- 
served in Faraday wave experiments, which result from 
the superposition of two frequencies. Although formu- 
lated primarily to describe those patterns, this model was 
also shown to describe other patterns, including regular 
square crystal lattices in two dimensions with the choice 
Qi/qo = V^, which couples (10) and (11) RLVs. 

The present introduction of the parameter ri in the 
form Q provides the additional flexibility to change the 
relative stability of different crystal structures. This is 
because in the limit fi 3> (/q' ^^is form reduces formally 
to the original Swift-Hohenberg form ([s]) after a simple 
rescaling of the parameters. Hence, as ri is increased the 
contribution of the second gi-mode becomes less signif- 
icant in comparison to the first gg-mode. Consequently, 
as ri is increased from zero, the crystal structure favored 
by the two-mode interaction becomes metastable with re- 
spect to the one-mode structure. This added capability 
to model the coexistence of two different crystal struc- 
tures, in addition to the coexistence of each structure 
with a liquid, should prove useful to model a wide range 
of phase transformations with a PFC approach. 

In the next section, we scale the parameters of the 
model to write the free-energy functional in a dimension- 
less form with only three parameters: e, which is the 
standard PFC model parameter analogous to tempera- 
ture that controls the size of the solid-liquid coexistence 
regions as a function of density, Qi = qi/qo, whose value 
is generally determined by the choice of crystal struc- 
ture, and i?i = ri/qQ controls the relative stability of 



the two-mode and one-mode structures (fee and bee, re- 
spectively). In this section, we also use a standard com- 
mon tangent construction to compute the phase-diagram 
in the plane of density and e for an illustrative choice 
of Ri = 0.05. The phase diagram exhibits regions of 
bcc-liquid and fcc-liquid coexistence for small and large 
epsilon, respectively. The size of the fcc-liquid coexis- 
tence region depends generally on Ri. For ri = where 
Eq. Q reduces to the free-energy density of Lifshitz and 
Petrich [18], the analog phase-diagram only exhibits fcc- 
liquid coexistence, so that a finite ri is necessary for the 
phase diagram to exhibit both bcc-liquid and fcc-liquid 
coexistence. We demonstrate the feasibility of the ap- 
proach with some simulations of polycrystalline growth 
and (111) twin growth. A numerical computation of the 
(111) twin boundary energy for parameters of Ni is given 
in an appendix. The ability to model twin growth is im- 
portant for solidification modeling since twins can dra- 
matically alter both eutectic pL9j 20: and dendritic pT] 
microstructures. 

In section HI, we carry out an amplitude expansion of 
the bulk free-energy density in the small e limit. This 
expansion exploits the property that, with the scaling 
i?i = eR, the amplitudes of the (111) and (200) den- 
sity waves scale as. As ^ Ae^^"^ and Bg ^ Be^'^, re- 
spectively, while the density difference between solid and 
liquid scales ~ e^/^. Therefore, this density difference 
can be neglected in the small e limit and the bulk free- 
energy density can be expressed solely in terms of those 
amplitudes. As required for solid-liquid coexistence, the 
free-energy density has minima in the {A,B) plane cor- 
responding to liquid {A — B = 0) and fee solid (finite A 
and B values that depend on R) . By comparing this form 
to the free-energy density for a single amplitude of bee 
density waves (corresponding to (110) RLVs), we identify 
different regions of relative fee and bcc stability, which 
explains the phase diagram computed in section II. 

In section IV, we discuss how to determine the two- 
mode PFC model parameters to relate them quantita- 
tively to different materials. We follow essentially the 
same approach developed by two of the authors for the 
standard PFC one-mode model for bcc ordering [T3] . For 
bcc, the parameters were completely determined by fit- 
ting three parameters: (i) the peak value of the liquid 
structure factor, S{qo), where go = |A'iio|, (ii) the sec- 
ond derivative of the fourier transform of the direct cor- 
relation function at this peak, C"{qo), and (iii) the solid 
density wave amplitude uhq. For fee, all the parame- 
ters except i?i are determined by the same fit, where 
Qn = 1-^111 1- (The shape of the structure factor at 
qi = |iir20o| is not realistically modeled given the lim- 
ited number of model parameters.) i?i then determines 
the ratio U20o/'"iii of the (111) and (200) solid ampli- 
tudes, which can be varied to alter the relative stability 
of fee and bcc. 

In section V, we derive analytical expressions for the 
three independent elastic constants of a cubic material, 
Cu, C12, and C44, for both the standard one-mode PFC 
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model ([3| and the present two-mode model Q. We 
use a brute force approach that consists of calculating 
to quadratic order the change of solid free-energy den- 
sity, modeled by a one- or two-mode approximation for 
bcc and fee, respectively, due to small dilation or shear 
transformations of the unit cell. We have checked that we 
obtain identical expressions to those derived recently by 
Spatschek and Karma for general lattices using an am- 
plitude equation framework |22j . which provides a non- 
trivial self-consistent test of our calculations. For the 
one-mode bcc model ([3]), the elastic constants are 
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Ci2 — C4 



noksT 



C"{qo)qWn„ (5) 



where go = l-f^iiol- For the two-mode fee model (|4 

C7n = -^^C"(go)g„^ (^?n + 4«Lo) , (6) 



and 



Ci2 = C44 = -^^^C"(<joKn, (7) 



where go = I ^111 1 Ri — for simplicity. 

Using values of C"(go) and density wave amplitudes 
from molecular dynamics simulations for parameters of 
bcc Fe and fee Ni, we find that the above expressions give 
reasonable estimates of elastic constants (e.g., Cn « 90 
GPa for one-mode bcc PFC model compared to Cn w 
128 GPa in MD Fe and Cn « 106 GPa for the two- 
mode foe PFC model compared to Cn « 155 GPa in MD 
Ni). The predicted values generally tend to be lower than 
the constants computed from MD simulations, but such 
discrepancies are to be expected given the PFC models 
are based on one or two modes. 

The analytical predictions for the elastic constants al- 
low us to draw two important general conclusions per- 
taining to the development of PFC models for different 
crystal structures and to the method used to determine 
the parameters of those models. 

The first conclusion, which follows directly from Eqs. 
^ and ([7]), is that the presence of the second mode, 
which corresponds to [200] density waves, is essential to 
obtain a physically meaningful set of elastic constants for 
fee. Without this second mode (u2oo = 0), Eqs. ^ and 
([7]), predict that Cn = C12 — C44. This implies that the 
tetragonal shear modulus C — (C12 — C22)/2 vanishes, 
and that the system is mechanically marginally stable. 
Of course, these analytical expressions for the elastic con- 
stants neglects the contributions of higher modes that are 
present in a full solution of the PFC equations. However, 
those higher modes are generally small for the small val- 
ues of e corresponding to Fe and Ni parameters. There- 
fore, the contributions of those modes will generally be 
small and will not change qualitatively this picture. 

While it is in principle possible to select energetically 
different crystal structures in the PFC model with the 
addition of other nonlinearities in the free-energy density 



(such as |V0|^ and V0p) [53], this approach will be of 
limited applicability for crystal structures like fee where 
one mode does not sufhce to produce the correct elastic 
properties. This is also true for simple cubic lattices and 
two-dimensional square lattices. The latter are briefly 
examined in section VI by coupling (10) and (11) density 
waves. 

The second conclusion, which is general, is that the 
elastic constants are uniquely determined once the phase- 
field model parameters have been fitted to the peak liquid 
structure properties, which fixes C"{qo), and the solid 
density wave amplitudes, as in the approach of Wu and 
Karma summarized above. This also fixes the value 
of the elastic bulk modulus 



K ^ 



Cn + 2Ci2 



(8) 



In general, the bulk modulus can also be defined from 
the thermodynamic relation 



K = V 



(fiF_ 



drfi 



(9) 



where F is the total free-energy, V is the volume, and 
n — N/V is the number density. The second equality 
in the last equation can in principle be used to compute 
the bulk modulus directly from the PFC solid free-energy 
curve {F/V versus n), without computing the elastic con- 
stants. For a perfect crystal without vacancy, Eqs. ([s]) 
and (9]) should in principle predict the same bulk modu- 
lus. However, the two definitions can give different pre- 
dictions for the PFC model because the number of atoms 
per peak of the crystal density field is not constrained to 
unity. While the average number of atoms per peak will 
also differ from unity in a real crystal with vacancies, 
thereby altering the open-system elastic constants |24) . 
the vacancy concentration is generally very small even at 
melting. How to meaningfully relate the predictions of 
Eqs. ([s]) and ^ for the bulk modulus is unclear in the 
PFC approach that, by construct, does not use a realistic 
description of the crystal density field, and also does not 
model vacancy formation explicitly. 

Despite these limitations of the PFC approach, Jaati- 
nen et al. |15) have recently proposed a modified one- 
mode PFC model to remedy the fact that, for the stan- 
dard one-mode PFC model with the free-energy density 
([3]), the bulk modulus predicted by Eq. ^ is several 
times smaller than the experimental value for parame- 
ters of bcc Fe. Their model yields a value of the bulk 
modulus computed through Eq. ([9]) that is in better 
agreement with experiment and also gives an improved 
prediction of the density difference between solid and liq- 
uid. It gives similar predictions of crystal-melt interfacial 
free-energies for bcc Fe as obtained previously by Wu and 
Karma using the standard one-mode model [H] . 

In the light of Eq. ([5]), it is apparent that any one- 
mode model that fits the correct peak structure factor 
properties and solid density wave amplitudes should pre- 
dict the same elastic constants. This is consistent with 
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the fact that Eq. ([s]) predicts a shear modulus C44 « 45 
GPa for the standard one-mode model of bcc Fe, which is 
reasonably close to the value C44 « 53 GPa estimated by 
Jaatinen et al. |15' from numerical shearing experiments 
in their model for similar input parameters. 

Since elastic constants are a major determinant of 
grain boundary energies and long-range interactions be- 
tween crystal defects, reproducing those constants, and 
hence the bulk modulus predicted by Eq. appears 
essential for modeling microstructural evolution. Also re- 
quiring that Eq. ^ predicts the correct bulk modulus 
using the solid free-energy curve may appear desirable. 
However, the motivation for doing so in the context of 
simple PEG models is somewhat less clear given the lack 
of realism of the crystal density field and the fact that 
Eqs. ([5])-([7| predict reasonable values of the elastic con- 
stants. In fact, any one- or two-mode model with the 
same peak liquid structure factor properties and density 
wave amplitudes will predict essentially the same elastic 
constants associated with the free-energy cost of lattice 
distortions, and also the same interfacial energies as can 
be inferred from amplitude equations [14 . Since those 
elastic constants and interfacial energies are the quanti- 
ties that matter most for modeling microstructural evo- 
lution in a PEC context, we have not found it necessary 
to formulate the two-mode PEG model in such a way 
that the bulk modulus is also correctly predicted from 
the solid free-energy curve using Eq. ([9|. Accordingly, 
we follow essentially the same approach outlined in Ref. 
[14j for determining the PEG model parameters. 



II. PHASE-FIELD CRYSTAL MODEL 
A. Basic equations and scalings 

The PEG equations have the standard form for con- 
served dynamics 



dr 



(10) 



where is the free-energy functional defined by Eq. ( 2 ) 
with the free-energy densities given by Eqs. ^ and (4) 
for the one- and two-mode models, respectively. To min- 
imize the number of parameters, it is useful to rewrite 
the equations in dimcnsionless form. Eor the two-mode 
model, we define the dimcnsionless parameters 
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Qi 



Mr 



ri 
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(11) 



(12) 



(13) 



where we set Qi = |J?2oo|/|^iii| = \/4/3 for fee (i.e. 
Qi equal to the ratio of the magnitudes of the (200) and 
(111) RLVs). We also define the dimensionless variables 
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F = 



J". 



(14) 



(15) 



(16) 



(17) 



Substituting the above definitions into Eqs. ^ and Q 
yields the dimensionless form 



'dt 



,6F_ 



(18) 



(19) 



with the free-energy functional 
and free-energy density 



/ = ^ H + (V^ + 1)^((V2 + + R,)] ^ + ^, 

where we have dropped the prime symbol on the dimcn- 
sionless spatial coordinate vector r' for brevity. Even 
though most of the paper focuses on the two-mode model, 
we also compute in section V the elastic constants for the 
standard one-mode PEC model. Eor this model, we use 
the same scaling as in Ref. [H] with the parameter 



Mr 



and dimcnsionless variables 



Mq 



t = FAggr, 



F 



(21) 



(22) 



(23) 



(24) 



where r' is defined by Eq. (|l4[ ). Substituting the above 
forms into Eqs. ^ and ([3| yields (after dropping the 
prime symbol on r') the dimensionless form of the one- 
mode PEC equations (IT8| and (figl) with 



/=tH+(V' + l)']^ + 



(25) 
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FIG. 1: Phase diagram of the two-mode PFC model for Ri = 
0.05 computed using two-mode and one-mode expansions of 
the crystal density field for fee and bcc, respectively. 
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FIG. 2: Phase diagram of the two-mode PFC model showing 
only the fcc-solidus and liquidus for the case iii = where 
fcc-liquid coexistence extends to vanishingly small e. 



B. Phase diagram 

The phase diagram of the two-mode PFC model is ob- 
tained by computing the free-energy density as a func- 
tion of the mean density ■0 in solid and liquid, denoted 
by /s (■)/'), and fi{ip), respectively, and then using a stan- 
dard common tangent construction to obtain equilibrium 
values of in solid (^ps) and liquid {ipi). 

Since the density is constant in the liquid, /; is ob- 
tained directly from Eq. ( 20 ) 
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Ri 



4 



(26) 



For small e, the solid free-energy density can be well ap- 
proximated by only considering the contribution of the 
(111) and (200) RLVs. Accordingly, the crystal density 
field is expanded in the form 



V'(r) 



E 

i?/ = (200> 



i?,=(iii) 
Tp + 8As cos qx cos qy cos qz 
+2Bs (cos 2qx + cos 2qy + cos 2qz) , (27) 

where we have used the fact that all density waves have 
the same amplitude in the crystal {\Ai\ = Ag and = 
Bg) and the magnitude of the principal RLVs are unity in 
our dimensionless units so {q = The parameters 

Ag and Bg are solved by substituting Eq. (27) into Eqs. 



19| and ( 20 1 and by minimizing the resulting free-energy 
and B, 



F with respect to A 
the solid free energy density 



/s(V's 



This minimization yields 



4 (-6 + 3^2)^2 + 3 + 3^2 ^ ^^^2 

45 

+72^sA^gBg + lUAlBl + 54A^ + —B^ 

+ + + (28) 

where Ag and Bg are themselves functions of ip. The co- 
existence densities V's and are computed numerically 
using the standard common tangent construction, which 
consists of equating the chemical potentials f'g{ips) = 
fliipi) = ^J.E and grand potentials /s(V's) - M-eV's = 
fi{ipi) — fJ-Eipi of the two phases. It is also necessary 
to compute the solid free-energy curve for bcc since the 
latter can have a lower free-energy than fee for some re- 
gions of the phase diagram. The bcc free-energy density 
was obtained by expanding the crystal density field us- 
ing a one-mode approximation, which only involves (110) 
RLVs as in Ref. [Mj , and substituting this expansion into 
the two- mode model defined by Eqs. (19) and (20). 



An example of the phase diagram for Ri — 0.05 is 
shown in Fig. [l] where we also show for completeness 
the hexagonal and stripe phases. As desired, we obtain a 
large e range of fcc-liquid coexistence. For small e, how- 
ever, bcc becomes favored over fee. A common tangent 
construction using fee and bcc free-energy curves shows 
that the density range of bcc-fcc coexistence is extremely 
narrow for small values of e and cannot be resolved on 
the scale of Fig. [T] As will be explained later in sec- 
tion [III C[ the range of e where bcc is favored depends 
on the value of In the limit Ri ^ 1, the two-mode 
model reduces to the standard one-mode model after a 
simple rescaling of parameters, which can be easily seen 
by comparing Eqs. (20) and (25). Hence, increasing Ri 



reduces the contribution of the second mode. Conversely, 
reducing i?i increases the contribution of this mode and 
tends to favor the fee structure, which extends to smaller 
e for smaller In the extreme case where i?i ~ 0, 
the region of fcc-liquid coexistence extends all the way to 
vanishingly small e as shown in Fig. [2] 





FIG. 4: Simulation of an equilibrium coherent (111) twin 
boundary for 7? = 0, e = 0.00823, and tp = -0.06269. 



C. Numerical examples 



We now demonstrate the feasibility of the model with 
some numerical examples of fee polycrystalline growth 
and (111) twin growth. The PFC conserved dynamics 
governed by Eq. (181 with the free-energy defined by 
Eqs. (19 1 and (20) was solved using the semi-implicit 



FIG. 3; Simulation of polycrystalline solidification starting 
from three seeded fee crystals in a supercooled liquid. The 
system is fully periodic, and the snapshots are taken at dimen- 
sionless times t = 5 x 10^,3 x 10^, and 10^. The parameters 
are i? = 0, e = .00823, and ip = -0.06. 



pseudo-spectral scheme given by Eq. (A2) in Appendix 
A of Ref . [25j . We used the parameters R = and e = 
0.00823 obtained from our fit of pure Ni presented later 
in section IV, together with the grid spacing Ax = Ay — 
Az = 27r'\/3/16, which determines the number of Fourier 
modes, and the time step At — 0.5. For this value of 
R and e, the computations presented in the next section 
show that the size of the solid-liquid coexistence region is 
extremely small, i.e. ips — ipi is two orders of magnitude 
smaller than {ips + ipi)/2 as can already be seen from the 
phase diagram in Fig. [2] and ^p^^ i^i ^ —0.0627. 

The first example in Fig. [3] shows the growth of small 
fee crystallites of different orientations for a value of 
-0 — —0.06 > -ips that is well inside the stable fcc-solid 
region of tlie phase diagram. The crystallites grow as ex- 
pected until they collide to form grain boundaries. The 
second example in Fig. |4] shows a (111) twin crystal for 
a value of ■0 = —0.06269 at coexistence and for a system 
size chosen such that a twin crystal with two stacking 
faults fits perfectly the periodic boundary conditions in 
all directions without any liquid present. A computa- 
tion of the excess free-energy of this twin boundary given 
in the appendix to this paper yields a value of approx- 
imately 30 mJ/m^ that falls within the range of values 
typically reported in the literature for fee metals. Fig. [5] 
then shows the growth of the same twin crystal in a super- 
cooled liquid for a much larger system with tp = —0.06. 
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FIG. 5: Simulation of the growth of a twin crystal in a super- 
cooled liquid for i? = 0, e = 0.00823, and i> = -0.06. 

III. AMPLITUDE EQUATIONS 
A. Scalings 

In this section, we analyze in more detail the proper- 
ties of the model by expanding the free-energy in terms 
of the amplitudes of density waves. In the one-mode bcc 
case analyzed in Ref. |14j . a similar expansion exploited 
the fact that the amplitude of (110) density waves scales 
as C""^/^ in the small e limit. In the present case, the ex- 
pansion is rendered more difficult by the presence of two 
different sets of density waves with amplitudes As and 
Bg corresponding to (111) and (200) RLVs, respectively. 
Therefore, it is not a priori obvious how As and Bg should 
scale in the small e limit. If Ri is kept constant, the bcc 
structure turns out to always be favored in the small e 
limit as apparent in the phase diagram of Fig. [T] Con- 
sequently, a small e amplitude expansion that captures 
the fee structure cannot be carried out at fixed i?i . How- 
ever, if i?i is decreased proportionally to e by imposing 
the additional scaling Ri = eR, both As and Bs scale as 
e^/^, thereby making a rigorous expansion possible. This 
expansion may seem artificial since the phase diagram of 
Fig. [l] is computed at fixed Ri. However, as we show 
below, the results of this expansion can be used to un- 
derstand the small e structure of the phase diagram, in 
particular the relative stability of fee and bcc. 

To demonstrate the feasibility of this expansion, we 
first analyze fee-liquid coexistence for small e with the 
scaling Ri — eR. The equilibrium densities are calcu- 
lated using the common tangent construction described 
in the previous section. To make the dependence of the 
coexistence densities on e explicit, we make a log-log plot 
of the mean coexistence density -0* = ^(V'/ + ips) versus e 




e 

FIG. 6: (Color Online) Plots of ijj* = |(V's + ipi) versus e 
for different values of R, with t/js and calculated from the 
common tangent construction. Fits to the numerical results 
of the form tp* = Vct^^^ yield for R ^ 0, ipc = -0.6901, for 
i? = 3, Vc = -0.6578, and for R = 5, ipc = -0.6396. The 
solid black line has a slope of exactly 1/2 on this log-log plot 
showing that the mean density ~ e^^^. 

for three different values of R. The results in Fig. [6] show 
that the mean coexistence density scales as e^/^. Next in 
Fig. [Tj we show a log-log plot of the density difference be- 
tween solid and liquid versus e for the same three values 
of R. The results show that ■0s — ^/i; ~ e^/^. Together, 
these two log- log plots show that, in the small e limit, the 
two-mode PFC model exhibits a weak first-order freezing 
transition where the size of the solid-liquid coexistence 
region is at the order of e'^/^ that is much smaller than 
the mean value of the density ^ e^/^. 



B. Free-energy functional 

The above scalings suggest that we can expand the 
crystal density field in powers of e^/^ as 

V'(f) = ,/>o(f)ei/2+^i(f)e + v^2(f)e3/2^... ^ (29) 
and expand accordingly the average densities 

0/ = V/oe'/' + V'/ie + V'/2e'/' + • • • , (30) 

and 

iPs = V'.oe'/' + 0.ie + V'.2e^/' + • • • , (31) 

in the liquid and solid, respectively. The numerically 
determined scaling relations (-0; -f -05)72 ^ e^^^ and ■0s — 
V'z ~ e'^^^ then imply that 

■0/0 = 0so = i>c, (32) 
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Since (V^ + 1)^(V^ + Qi)^V'2 gives a vanishing contri- 
bution for all density waves associated with sets {Ki} 



iKi 



and 



must 



and {K'j\^ all remaining terms 
balance each other in order for a solution of Eq. (|39|) to 



exist. For example, the condition that the coefficients of 



balance each other yields 



(-1 + 3Ve')^?n + (3|^?„p + 6|^;ij2 + 6KiP 
+6|A?„|2 + 6|i?°ool' + 6|So°2ol' + 6|i?o°o2l')^ni 
+6i/'c(^lii-B2oo + ^111^020 + ^111^002) 

+""^111"^111^111 + O^lll-°200-°020 

+6^11-^200-^002 + 6^111-8002-^020 = 0' (40) 



FIG. 7: (Color Online) Plots of tps — tpi versus e for the same 
values of R as in Fig. [sj with ips and calculated from 
the common tangent construction. The solid black line has 
a slope of exactly 3/2 on this log-log plot showing that the 
density difference between solid and liquid ~ e^^^ . 



and 



(33) 



Next, to carry out the amplitude expansion, we start 
from the equilibrium equation 5F /5ip = /zg, where /i^; is 
the equilibrium value of the chemical potential. With F 



defined by Eqs. (19) and J20p, we obtain 



/iij = -eijj + (V^ + 1) [(V^ + Qlf + ei?] V + V-^- (34) 
We substitute the small e expansion of the density field 



( 29 1 into Eq. ( 34 1 and collect terms with the same power 
of e. We find at the order e^/"^ 



which has the solution 



iKj -r 



(35) 



(36) 



where the summations are over (111) and (200) RLVs, 

respectively, and \Ki\ = 1, | = ^4/3, in our scaled 
units. At order e, we obtain 

(v^-f i)2(v2 + g2)2^i = 0, (37) 

which has the solution 

^1 = ^ Ale^'^-^^ + ^ B]e^'^^' \ (38) 

i 3 

and collecting the terms at order e'^/^ yields 

-Vo + (V^ + 1)2(V2 + Q2)2^2 ^ ^(v2 ^ ^)2^^^ ^ (^^)3 

= -V'c + Q\i^i2 + mc + (V-c)^. (39) 



and requiring that the coefficients of e'^^"" '' balance each 
other yields in turn 



(-1 + 3V',2 + R[-Ql + lf)Bl^^ + (6|A;nP + elAjjJ 



0201 



-f6|A;,lp-f6|A;„P + 3|i?°ooP + 6|i3i 
+6|i?o"o2nS2^o + Hc{Al,-,Al^, + ^?ii^?ii) 

+"l-D020^111^111 + -°002^111^111 

,nO 4O 4O , nO 4O 4O \ n 
+ -°020"^111"^111 + -°002^111^11lJ — 



(41) 



The above solvability condition must be satisfied inde- 
pendently for each reciprocal lattice vector. This yields 
a set of fourteen coupled amplitude equations that are 
straightforward to obtain. From those amplitude equa- 
tions, it is useful to express the free-energy of the system 
measured from its constant value in the liquid, defined 
as the difference I^F^^ , as a functional of the density 
wave amplitudes and S?. This quantity can be ex- 
pressed solely in terms of the amplitudes of density waves 
owing to the property that the size of the coexistence 
region (~ e^/'^) is much smaller than the mean density 
(~ e^/2) in the small e limit. Since the amplitudes are not 
conserved order parameters, the equilibrium state simply 
corresponds to a minimum of this free-energy. Hence the 
amplitude equations must be equivalent to 



5AF 



AE 



5A^* 



(42) 



and 



5AF 



AE 



SB?* 



^ 0. 



(43) 
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For the case where the amphtudes are spatially uniform, 
we obtain the free-energy density 



(-l + 37/^,^)(|A'/ 



111 



Kil' + I^ml 



(44) 



1-^200 I 



\A 



l^ml' + l^ml 



+6i\AU'\AU 
+\AW\A'inr 
+\A",,,mnr 



' + |SoVr) 
+ \AU'\AU' 
\AU'\AU' 



I -^020 I I -^002 



0201 
2\ 



MAini'mu 

+6MAUAlriBloo 

j_aO _ aO dO_ 
"r^lll^lll^020 



1^2001 

6(l^?ii 

- I^020P ■ 



2|r0 |2 
1^002! 
2 



I-S002P 



aO aO rO 
^lll^lll-°020 



_ 
002 



^^111^111^002 
T^^lll^lll-°002 
^^^111^111^200 
+"^lll^lll^lll^lll 

_i_fi/io_ AO__nO oo 

"r"^lll^lll-°200^020 
_lK/10__40 nO rO 
"T "^lll^lll-°002^020 

+6^111^111-^^00^020 

+6^111^111-^200^020 

+6^111^111-^200-^002 

+6^111^111-^002-^020 



A" - A"- Ry 

^111^111-^200 
/lO aQ tjO 
^111^111-^002 

A^- /10__rO 
^lll"^lll-°020 

aO aQ nO 
^111^111-^020 

aO _ aO_ rO ^, 
^111^111-^200 



fi/io 4" 40 40 
"^111^111^111^111 

ftAO /lO nO rO 
"^lll^lll-°200-°002 

rO 
002^020 



a aO_ aO rO rO 

""^lll^lll-°200-°002 

GA^n^^iii?; 

6^111^1Il-^20n^i 

6^;ii4iii?r 



rO 
200-^020 

dO 



rO- 
002^020- 



In general, the above free-energy density is a multi- 
variate function of the fourteen amplitudes A° and i?f of 
different density waves, and cannot be represented graph- 
ically in a simple way. However, the free-energy barrier 
between solid and liquid can be made explicit by assum- 
ing that all the (111) and all the (200) density waves have 
the same amplitude (i.e., A'^ = A and = B, respec- 
tively). In this isotropic approximation (see jH] for the 
bcc analog), the free-energy density becomes 



A/^^ = 4(- 



-3(-l 
^54^4 



45 



1 

2 

-B" 



- 3^Ac')A2 
R{-Ql + 

+- UAA'^B^ 



0.15 



l^?iil 




0.05 



f)B^ 



72il:^A^B. (45) 



This expression can also be obtained by evaluating di- 
rectly the difference between the solid and liquid free- 



energy densities 
given by_Eqs 
tutions tps = Tpi 



A/ 



AE 



e {fs- fi), with fi and 



(|26[ ) and ( 28 ) , respectively, and the substi- 
ipce^^A, = and B, = Be^'"^. 

For the parameters i? = and tAc calculated from the 



FIG. 



(451 



8: (Color online) Free-energy landscape defined by Eq. 
as a function of the amplitudes A and B of (111) and 
(200) density waves, respectively, for R = and the corre- 
sponding coexistence value V'c = —0.6901 where the solid and 
liquid minima have the same height. 



common tangent construction, we plot in Fig. [8] the free- 
energy landscape as a function of amplitudes A and B. 
The free-energy landscape exhibits two minima that cor- 
respond to the stable liquid and solid phases. The above 
amplitude equation calculation shows that the two-mode 
PFC model describes well solid-liquid coexistence with a 
well-defined free-energy barrier between solid and liquid. 

We have only treated here the case where the ampli- 
tudes are spatially uniform to characterize the bulk free- 
energy landscape. A more general free-energy functional 
that includes gradient terms would be necessary to treat 
the case where the amplitudes are spatially varying. Such 
a functional could then be used to compute the excess 
free-energy of the solid- liquid interface and its anisotropy, 
as done previously for bcc jl4j. Those computations will 
be presented elsewhere. 



C. Relative stability of fee and bee 

So far, we have only examined the possibility of fcc- 
liquid coexistence. However, the phase-diagram of Fig. 
[T] shows that bcc can have a lower free-energy than fee for 
small enough e if i?i is finite. We now use the amplitude 
equations to understand the relative stability of fee and 
bcc. As a first step, it is useful to re-examine the scaling 
of the mean density that is controlled by the parameter 
ipc- We computed previously the equilibrium solid and 
liquid densities using the common tangent construction, 
from which we obtained the scalings {'ips+'>Pi)/'2 = i/'cC^^^, 
which defines ipc, and ips ~ '4'i e"^/^, which shows that 
the size of the density difference between solid and liq- 
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Stable Fee, Metastable Bee 



Stable Fee, Unstable Bee 



R 



FIG. 9: (Color Online) i/jc as a function of R showing different 
ranges of bcc or fee metastability, stability, or instability. 



uid can be neglected in the small e limit. We can also 
compute ipc more directly from Eq. (451 by requiring 



dAf 



AE 



dAf 



AE 



OA 



dB 



0, 



and 



A/^^ = 0, 



(46) 



(47) 



with all the above relations evaluated at the equilibrium 
values of A and B in the solid. Eq. (46) stems from the 



requirement that the solid amplitudes must correspond 
to a free-energy minimum, which fixes those amplitudes 
uniquely as functions of ^pc■ Eq. (471, in turn, is the re- 



quirement that the free-energies of solid and liquid must 
be equal in equilibrium, which fixes ijjc uniquely for a 
given R. A plot of "fAc versus R obtained in this way 
using Eqs. (|46| and ^ is shown in Fig. |9] 

The relationship between ipc and R, denoted by ^c{R), 
can now be used to assess the relative stability of bcc and 
fee. To obtain an analogous expression to Eq. (45) for 



bcc, we substitute into the two-mode free-energy func- 
tional defined by Eqs. (19) and (20), the one-mode ex- 



pansion of the bcc crystal density field in terms of prin- 
cipal set of (110) density waves 



■4>{r) ~ £^^^"00 + 4e^^^A(cos qx cos qy 



cos qx cos qz + cos qy cos qz) 



(48) 



where q = l/\/2. We obtain 



Af^^^ = 6(3^/;," - l)A^ + 48V'c^^ + 135^*, (49) 

where A now denotes the amplitude of (110) density 
waves and we have used the subscript "bcc" to distin- 
guish this free-energy difference between bcc and liquid 



li3 o 
<l 



0.06 



0.04 



0.02 




-0.02. 



FIG. 10: Plots of free-energy density of bcc relative to the 
liquid, A/i^f, as a function of the amplitude A of (110) den- 
sity waves for different values of R. When the value of A/^f 
corresponding to the solid bcc free-energy minimum is nega- 
tive, bcc is favored over fee. As R decreases, bcc first becomes 
metastable with respect to fee and then unstable as the local 
solid free-energy minimum disappears. 



from the one between fee and liquid, A/"*^, defined by 
Eq. ( 45 ) . By definition, Af^^ = for solid fee in equi- 
Therefore, 



librium with the liquid. Therefore, to assess the relative 
stability of bcc and fee, we can plot Af^f defined by 
Eq. (49) versus A and check if the value correspond- 



ing to the solid bcc minimum is above (below) zero in 
which case fee (bcc) has a lower free-energy than bcc 
(fee). Such plots shown in Fig. 10 show that bcc be- 



comes metastable with respect to fee and then unstable 
(with the disappearance of the local solid bcc free-energy 
minimum) as R is decreased. A detailed study as a func- 
tion of R shows that bee first becomes metastable for 
R < Rc where Rc ~ 2.68 and then unstable as R is de- 
creased below a second threshold value (~ 1-43), giving 
rise to the three different stability regimes as a function 
of R shown in Fig. [9j Translated in terms of the phase 
diagram constructed at fixed this implies that bcc 
becomes favored over fee when e < Cc where 



Rc 



(50) 



For i?i = 0.05, the above expression predicts ec ~ 0.019 
that is in good quantitative agreement with the phase 
diagram of Fig. [T] As i?i increases, ec increases and 
the switch from stable bee-liquid to fee-liquid coexistence 
moves to higher values of e in the phase diagram. 



IV. PARAMETER DETERMINATION 

In this section, we derive expressions to relate the two- 
mode PFC model parameters to material parameters by 
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extending our previous approach for bee |14j . As a first 
step, we match the peak Uquid structure factor properties 
of the two-mode PFC model to the standard expression 
from classical DFT. The expression for the PFC liquid 
structure factor is obtained by varying -0 around its liquid 
value, ij} — ijji + Sijj, and evaluating the corresponding 
variation AJ- of the dimensional free-energy difference 
between solid and liquid using Eqs. ^ and Q, and the 
relation (15 1 between (p and iJj. Dropping terms of S^p 



higher than quadratic order, we obtain 



PFC 



Ml 



dr 



Sip 



(51) 



Substituting the Fourier transform, 

dk 



Si: ^ 



(27r)3/^ 



(52) 



we obtain 



Af/Q [ f dkdk' StpkSijjk 
^ ~Y j J ^ 2~ ^ 

+ X{-k^ + qlf{{-k^ + qlf + r, 



a + 3ijfXql 



^ g J * 2 ' ° 

+ X{^k^ + qini-k' + qir + n)]. (53) 

A second expression for the free-energy of a spatially in- 
homogeneous liquid is obtained from classic DFT 



knT 



DFT 



Sn{r) 



S{r-r') 



no 



drdr 



Ci\f-r^\) 



Snir'), (54) 



where 



and 



Sn{r) — n{r) — Uq — 5(j){r) 



C(k) ^ no / drCi\r\)e 




5ip{r), 



-ik-r 



(55) 



(56) 



is the Fourier transform of the direct correlation function. 
Fourier transforming again, we obtain 

A^DFT ^ M /■ dkS^JkSi^., [1 - C{k)] . (57) 
g 2no J 



Equating A J"pfc = AJdft and using the expression for 
the liquid structure factor S{k) = 1/(1 — C(fc)), we obtain 



Sik) 



knT 



(58) 

By evaluating the above expression at the peak of the 
liquid structure factor, we obtain 



a + SXq^i^f = 



ksT 
noS{qo) ' 



(59) 



or, using Eq. (11) and the relationship tpi — ?/'ce^^^, 

' no5(go)Ago'(l-3^2)- 
A second relation is now needed to determine e and A 



(60) 



independently. To obtain it, we substitute Eq. (58) into 



the relation C{k) — {S{k) — 1)/S{k) and compute the 
second derivative of C{k) evaluated at the peak of the 
liquid structure factor to obtain 



A = - 



kBTC"{qo) 



^noqlil 



eR) 



Eqs. (60 1 and (61) combined now give 



9(g25(go)C"(qo)(l-3V'?)-8i?)' 



and 



~9kBTC"iqo) 



QkeTR 



no5(9o)(l-3V'c')'Z§' 



(61) 



(62) 



(63) 



In addition, the relation ([55| between the real and di- 
mensionless densities expresses 



Aq|^ 

9 9 
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(64) 



in terms of the solid amplitude of the first go-mode. 
The two solid amplitudes Ag and Bg can be computed for 
a given R by using the scaling relations As — e^^'^A and 
= e^/'^B where A and B are the equilibrium values of 
the scaled amplitudes in solid. The latter are obtained, 
together with ipc (Fig. [9]), by using the conditions (46) 
and ^ with Ai^^^ defined by Eq. 

For a given R, Eqs. (|62|, (|63]), and (^ fix the three 
parameters e. A, and g of the PFC model uniquely in 
terms of peak liquid structure factor properties, S{qo) 
and C"(go)j where qo = |i?iii| here, and the solid density 
wave amplitude Us — itm. This still leaves the freedom 
to vary R within the range where fee is stable with re- 
spect to bee (Fig. [9|. Varying R changes the shape of 
the liquid structure factor as shown in Fig. [12] and de- 
creasing R below some threshold produces a second peak 
at qi = |i^20o|) and generally increases the contribution 
of the second mode. Thus decreasing R increases the am- 
plitude of the second mode U200 as shown in Fig. 11 For 
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FIG. 11: (Color online ) P lots showing the variation with R 
of e calculated by Eq. (62 1 and the ratio Bs/As of the solid 



amplitudes of the se con d and first modes calculated by the 
conditions (lief and (|47| with AF^^ defined by Eq. ((45 1. 



simplicity, we used the value R = that yields a reason- 
able fit of this amplitude for pure Ni. The other input pa- 
rameters computed by Hoyt [37] using the EAM potential 
of Foiles, Baskes and Daw [28] (FBD) are given in Table|T] 
The density wave amplitudes are calculated using the re- 
lation Ui = ex-p{— Kf / Aa) , which assumes that the crystal 
density field is a sum of Gaussians centered around each 
fee lattice site. The value of a is obtained from the ex- 
pressions for the root- mean-square displacement of atoms 
in the solid a/ < |r]^ > = 3/ (2a) derived from this den- 
sity field. For the value ^< |f|2 > - 0.298 A from MD 
simulations, we obtain um — exp{—Kfii/4:a) = 0.6639 
and U200 = exp{- Kf^^^ / Aa) = 0.5791. 

It should be noted that, with the present fitting pro- 
cedure, the two-mode PFC model only reproduces the 
correct shape of the main peak of the liquid structure 
factor. The second peak is spurious and is only used 
to increase the amplitude of the second mode to some 
desired value. Since the second mode is critical to ob- 
tain solid-liquid coexistence, the lack of realism of the 
structure factor outside of the first peak is a limitation 
of the present two-mode model. The liquid structure 
factor could in principle be made more realistic by shift- 
ing the second peak to larger wavector and reducing its 
amplitude, which would couple the principal (111) RLVs 
to other sets such as (222) and (311). However, larger 
fc-modes require a finer mesh and are computationally 
more costly to resolve. Whether such a fit would offer 
specific advantages remains to be investigated. 



- MD Nickel Structure Factor 

-R=0 

-R=1 

. R =2.6 




FIG. 12: (Color online) Liquid structure factor of the PFC 
model and from MD simulations of pure Ni |27j . 



TABLE I: Input parameters for the PFC model computed 
from MD simulations of pure Ni [27] using the FBD EAM 
potential [28] and corresponding PFC parameters. 



MD input parameters 


Value 


PFC parameters 


Value 


go = Kill (A-i) 


3.0376 


A (eVA") 


0.026 


no (A-') 


0.0801 


g (eVA^) 


8.53 


C"{qo) {A') 


-9.1579 




-0.6901 


S{qo) 


2.9898 


R 





Tm (K) 


1811 


e 


0.00823 


itiii 


0.6639 


itiii 


0.6639 


U200 


0.5791 


^200 


0.5136 



pare the results to MD computations of elastic constants 
at the melting point for parameters of fee Ni. For com- 
pleteness, we also carry out the same comparison for the 
standard one-mode PFC model for parameters of bee Fe. 
Following the same approach as in Ref. [5], we obtain 
the elastic constants by deforming the lattice from its 
ideal structure and computing the corresponding change 
of free-energy density. We consider three different defor- 
mations 



V. ELASTIC CONSTANTS 

In this section, we derive analytical expressions for the 
elastic constants of the two-mode PFC model. We com- 



= ^ix/{i + 0,y/a + 0,^/a 

iP^{f) = ij{x + £,y,y,z), 



-0), 



(65) 
(66) 
(67) 
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of the two-mode crystal density field We compute the 
change of free-energy density 



(68) 



where fs is the free-energy density of the unperturbed 
solid and Fi is the free-energy integrated over the per- 
turbed unit cell of volume Vj with 



I M^+i) ra(i+i) Mi+a 

V^ " W+WJo Jo L 

1 t-a /-all-?) 1-0.(1+^) 

% ^ v{i - e) Jo Jo 



^2 
^3 



1 

VJo 



a pa pa—^y 



J-iy 



where dV = dxdydz, a is the lattice spacing, and V = 
is the unperturbed unit cell volume. 



fee elastie eonstants for the two-mode model 



Using Eqs. (65)-(68l with the two-mode crystal den- 



sity field ^p{r) defined by Eq. (27) and the free-energy 
density f{fp{r)) defined by Eq. (20), we obtain the di- 
mensionless elastic constants 



2^11 



A/2 = (Cn - Ci2)e 



ia + l3)e, 
2 



A/3 = 



(a/9 + S)^^ 



where we have defined 



^(284 - 315Q? + 8igt + Bli?!)^^ 



(69) 
(70) 
(71) 

(72) 
(73) 

(74) 



We can set i?i w in the above expression since Ri <C 
1 for typical model parameters where the fee lattice is 
favored. With only the principle (111) RLVs (Bg = 0), 
all three elastic constants are equal Cn = C12 = C44 = 
2a/9, which gives a vanishing tetragonal shear modulus 
C' = {Cii - C'22)/2. The inclusion of the (200) RLVs, 
however, raises the value of Cn, which becomes 

.Br 



while leaving the values of C12 and 6*44 unchanged, 
thereby making C" finite as desired. Finally, convert- 
ing back to dimensional units using the relation Cij — 
[X^q^^ / g)Cij, we obtain 



Cn = --nQkBTC"{Kni)qluln ( 1 



-^200 



''111 



(75) 



TABLE IL Comparison of elastic constants at the melting 
point predicted by the two-mode PFC model for fee Ni and 
the one- mode PFC model for bcc Fe and MD simulations 1301. 



Quantity 


PFC fee 


MD fee 


PFC bcc 


MD bee 


Cn (GPa) 


106.6 


155.4 


90.0 


128.0 


C12 (GPa) 


31.4 


124.7 


45.0 


103.4 


C44 (GPa) 


31.4 


66.0 


45.0 


63.9 


Bulk Modulus (GPa) 
(Cn -f2Ci2)/3 


56.5 


134.9 


60.0 


111.6 



and 



C12 — C4 



- no fc_B TC" (Xn 1 ) <7o 1 1 • 



(76) 



The elastic constants computed with the parameters of 
Table |I] are compared to the predictions of MD sim- 
ulations in Table HIl The MD simulations for fee Fe 
and bcc Ni were carried out using the EAM potentials 
from Mendelev, Han, Srolovitz, Ackland, Sun and Asta 
MH(SA)2 [55]^ and Foiles, Baskes and Daw [55], respec- 
tively. The same EAM potentials were used to compute 
the input parameters for the PFC model and the elastic 
constants. The input parameters for Fe are the same as 
in Ref. [26] . The input parameters for Ni were computed 
by Hoyt [57]. The elastic constants for both Fe and Ni 
were computed by Foiles [30]. Their values at the melt- 
ing point are smaller than at zero temperature as shown 
by Foiles for a different Ni EAM potential [31 . 



B. bee elastic constants for the one- mode model 



For the one-mode model, we use again Eqs. (65l-(68) 
with the one-mode bcc crystal density field 

~ '(/'-f 4^s(cos qx cos qy+cos qx cos qz+cos qy cos qz). 

(77) 

where q — 1/^/2 and the free-energy density fiipir)) de- 
fined by Eq. ([25|). We obtain 



A/i 
A/2 

A/3 



2^" 



(Cn- 

CiA ^2 

2 ^ 



+ 3Ci2 ) i 



(^bcc /-I 



(78) 
(79) 
(80) 



where abac = 24A'^. This yields the diniensionless elas- 
tic constants Cn = 2C\2 = 2(744 = 8^s- Finally, 
converting back to dimensional units using the relation 
Cij ~ (y? I g)C ij ^ we obtain 

Cn = 2 C12 = 2 C44 = -nofcBTC"(ifno)'?o«?io- (81) 

The elastic constants computed with the input parame- 
ters of Table |III| for bcc Fe are compared to the predic- 
tions of MD simulations in Table |lll 
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TABLE III: Input parameters for the PFC model computed 
from MD simulations of pure Fe [26 using the EAM potential 
from MH(SA)^ 29 and corresponding PFC parameters. 



MD input parameters 


Value 


PFC parameters 


Value 


go = ifiio (A-i) 


2.985 


A (eVA^) 


0.291 


no (A-') 


0.0765 


g (eVA«) 


9.705 


C'iqo) (A') 


-10.40 


e 


0.0923 


S{qo) 


3.012 


Who 


0.72 


Tm (K) 


1771 






UllQ 


0.72 







VI. TWO-DIMENSIONAL SQUARE LATTICE 

As an example of application of the two-mode model 
to other lattice structures, we briefly examine the exam- 
ple of two-dimensional square lattices. Those lattices are 
obtained by coupling (10) and (11) density waves with 
Qi = as demonstrated previously by Lifshitz and 
Petrich [TS] for a modifed Swift-Hohenberg model that 
corresponds to the i?i = limit of the present two-mode 
model. The liquid free-energy density is given by 



/f^(V^0 = "(e-4-i?i) 



(82) 



and the solid free-energy density is obtained by substi- 
tuting the two-mode crystal density field 

ip{f) ~ -0 -I- 2 As (cos X + cos y) + ABg (cos x cos y) (83) 



into the free-energy functional defined by Eqs. (191 and 
(Eol with Qi = V2 , which yields 



/|«(V;.) = 2(-e + 3^f)A,2 + 2(-e + 3V;2^i?i)i?,2 



Ri 



1 



(84) 



For Ri = 0, we obtain ipc — —0.6782 numerically from 
a log-log plot of the mean equilibrium density versus e 
similar to Fig. [6] which is determined from the common 
tangent construction. The feasibility of the two-model to 
model polycrystalline solidification and grain boundaries 
is illustrated in Fig. [13] As for fee, the second mode 
turns out to be essential to obtain physically meaningful 
elastic constants. Following the same procedure as for 
fee in the last section (with deformations of the unit cell 
now constrained to the x — y plane) we obtain 



(cii + C12) e 

t2 



■/3s«)r 

t2 



(Cll-Cl2)C' = (a,,+A,,)e^ (86) 
2 ^ 



(85) 
(86) 

(87) 




FIG. 13: (Color online) Example of polycrystalline solidifica- 
tion for two-dimensional square lattices. The snapshots are 
at dimensionless times t = 10, 100, and 1000. The parameters 
are e = 0.15, Ri = 0, Qi = \/2, and i> = -0.23. 
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with 



asq = ^{l-2Q\ + Q\ + Ri)Al, 
I3sq = 8{58~37Ql + 5Qt + 5Ri)Bl (89) 



Asq = 8RiBg, 



Jsq 



2 (90) 
4(32-2ig2 + 3g^ + 3i?i)s2. (91) 



Again we look in the limit _Ri w 0. This yields the di- 
mensionless elastic constants 



Cll ^ asq + /3sq/2, 

(744 = 2(5. 



(92) 
(93) 
(94) 



These relations show that the one-mode crystal density 
field consisting only of a superposition of (10) density 
waves {Bs — 0) yield vanishing shear moduli, which be- 
come finite with the inclusion of the second mode. 



VII. CONCLUDING REMARKS 

In summary, we have presented a two-mode PFC 
model with a phase-diagram that includes different tem- 
perature ranges for bcc-liquid and fcc-liquid coexistence. 
The relative sizes of these ranges can be changed by vary- 
ing one model parameter that controls the relative mag- 
nitudes of the amplitudes of the two modes, correspond- 
ing to [111] and [200] density waves, respectively. We 
have shown that the free-energy landscape for fcc-liquid 
coexistence has a double- well structure with a finite free- 
energy barrier between solid and liquid in the plane of the 
amplitudes of the two modes. We have demonstrated the 
feasibility of the model with some numerical examples of 
fee polycrystalline growth and twin growth, as well as for 
two-dimensional square lattices. 

At a more quantitative level, we have determined the 
model parameters by fitting the peak liquid structure 
factor properties (S'((?o) and C"(qo)) and solid-density 
wave amplitudes as an extension of our previous study 
of bcc Fe [14] . Furthermore, we have derived analytical 
expressions for the elastic constants. With input values 
for those parameters from MD simulations of pure Ni, 
we have found that the PFC model elastic constants are 
in reasonable agreement with MD results given the sim- 
plicity of the model, which neglects the contributions of 
many other modes that are present in a realistic descrip- 
tion of the crystal density field. Those expressions also 
stress the necessity of having at least two distinct modes 
to obtain physically meaningful values of the elastic con- 
stants for fee in the physically relevant small e-limit of 
the PFC model, which is also true for square lattices. 



We have found that the standard one-mode PFC model 
also predicts reasonable values of the elastic constants 
for pure bcc Fe, and we have argued that any one- or 
two-mode model will predict similar elastic constants for 
bcc and fee with the same peak liquid structure factor 
properties and solid density wave amplitudes 

Finally, while the numerical examples focused on crys- 
tal growth, it might also be possible to use the two-mode 
PFC model to study the bcc/fcc martensitic transforma- 
tion, which has been modeled by other phase-field ap- 
proaches that make use of structural order parameters 
(see Refs. [351 [33] and references therein). The abil- 
ity to vary the relative stability of fee and bcc crystal 
structures, which was demonstrated here, should prove 
particularly useful for this application. 
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Appendix A: Twin boundary energy 

We computed the coherent (111) twin boundary en- 
ergy using the method put forth in Ref. [23], which ex- 
ploits the dependence of the free-energy on system size. 
We performed simulations for four different lengths Lz 
along the axis perpendicular to the boundary. By plot- 
ting the bulk free energy density / against the inverse 
of this length (Fig. 14 1, we then extracted the boundary 



energy from the slope of this plot using the relation 



It-u 



(Al) 



where /s(V') is the free energy density of a perfect crystal. 
This method turns out to be more accurate than com- 
puting directly the excess free-energy of the boundary for 
a fixed system size [25]. This calculation gives twice the 
boundary energy since there are two boundaries in our 
periodic system. We convert the result to dimensional 
units through, "/twin = {>^^<lo' /g)ltwin, which yields the 
reasonable value of 29.9 mJ/m^ for the same Ni param- 
eters used to compute the elastic constants. 
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